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In this work we study the effect of metachronal waves on the flow created by magnetically- 
driven plate-like artificial cilia in microchannels using numerical simulations. The simula- 
tions are performed using a coupled magneto-mechanical solid-fluid computational model 
that captures the physical interactions between the fluid flow, ciliary deformation and 
applied magnetic field. When a rotating magnetic field is applied to supcr-paramagnctic 
artificial cilia, they mimic the asymmetric motion of natural cilia, consisting of an ef- 
fective and recovery stroke. When a phasc-diffcrcncc is prescribed between neighbouring 
cilia, metachronal waves develop. Due to the discrete nature of the cilia, the metachronal 
waves change direction when the phase difference becomes sufficiently large, resulting in 
antiplectic as well as symplectic metachrony. We show that the fluid flow created by the 
artificial cilia is significantly enhanced in the presence of metachronal waves and that 
the fluid flow becomes unidirectional. Antiplectic metachrony is observed to lead to a 

considerable enhancement in flow compared to symplectic metachrony, when the cilia 
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spacing is small. Obstruction of flow in the direction of the effective stroke for the case 

of symplectic metachrony was found to be the key mechanism that governs this effect. 



1. Introduction 



The control of fluid flow in channels of micron-scale dimensions is essential for proper 
functioning of any lab-on-a-chip device. The fluid transport in microchanncls is of- 
ten pe rformed by downscaling conventional methods such as syring e pumps, microp- 



umps (jLaser fc Santiago 



2004; 



Jeon et al 



2000; 



Schilling et al 



20021 ) or by exploiting 



electro- magnetic fluid manipulation principles, as in electro-osm otic (jChen et al. 



Zeng et al 



2003 



20021 ) devices. In search for 



20021 ) and magneto-hydrodynamic (jWest et al. 
novel ways to propel fluids at micron scales, we let nature be our guide. Nature uses 
hair-like structures, called cilia, attached to the surfaces of micro-organisms, to pro- 
pel fluids at small length scales. The typical length of a cilium is 10 microns. Cilia 
beat in a whip-like asymmetric manner consisting of an effective stroke and a recovery 
stroke. Moreover, when m any cilia operate to gether, hydrodynamic interactions cause 



them to beat out-of-phase (IGueron et al. 



waves, and an enhanced fluid flow 



1997), leading to the formation of metachronal 



Satir fc Sleigh 



1990I ). The specific metachrony is 



termed symplectic (or antiplectic) when the metachronal wave is in the same (or op 
posite) direction as the effective stroke. The cilia on a Paramecium ex hibit antiplectic 



meta chrony, whereas the cilia on Opalina exhibit symplectic metachrony ([Blake fc Sleigh 



19741) . The asymmetric motion of natural cilia is due to the intricate interaction between 
the cilia micro-structure (axoneme) and the internal driving force generated by ATP- 
enablcd conformational changes of the motor protein dynein. It is a challenging task to 



design the artificial counterpart of natural cilia, by using external force fields for actu- 
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ation in order to mimic the asymmetric motion of natural cilia. An early attempt to 

create art ificial cilia was based on electrostatic actuation of arrays of plate-like artifi- 



cial cilia ( den Toonder et 



2008). Although effective flow and mixing were achieved, 
movement of these artificial cilia was not asymmetric as in the case of natural cilia. It 
was predicted using numerical simulations that an array of identical super-paramagnetic 
or permanently magnetic two-dimensional plate-like cilia can mimic the planar asym 



metr ic motion of natural cilia when exposed to a uniform magnetic field (jKhaderi et al 



2009f ). These magnetic plate- like cilia can be realised, for instance, by using polymer 



film s with embedded super-paramagnetic (or permanently m agnetic) nano-particles (see 



e.g. 



Fahrni et al 



2009; 



Belardi et al 



2010 



Schorr et al 



2010). In contrast with the plate- 



like cilia, rod-like structures that mimic the three-dimensional motion of nodal cilia to 



create fluid propu lsion have also been fa bricated (jVilfan et al 



Evans et al 



2003). In 



Sing et al 



2010: 



Shields et al 



2010 



(|20I0l ). a novel method of fluid propulsion based on 



magnetic walkers was prese nted. Artificial cilia bas ed on photo-actuation have also been 



realised in the recent past (jvan Oosten et al 



20091 ) 



In previous numerical studies we focused on the flow created by an array of s ynchronously 



beati ng plate-like cilia whose motio n is planar and asymm etric, in the absence (jKhaderi et al 



2009) and presence of fluid inertia (jKhaderi et al 



20101 ). It was reported that a substan- 



tial but fluctuating flow is created in the former, while in the latter the flow increases 
significantly as the Reynolds number is increased. In addition, the fluid flow can become 
unidirectional in the presence of fluid inertia. In this work we explore another aspect 
of natural ciliary propulsion using numerical simulations - the metachronal motion of 
cilia, by allowing the asymmetrically-beating artificial cilia to move out-of-phase. The 
out-of-phase motion of the cilia is achieved by applying a magnetic field that has a phase 
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lag between adjacent cilia. The existing literature on the metachronal motion of natural 

cilia could provide insights on the flow generated in the presence of metachronal waves. 

In the case of natural cilia the metachronal motion is analysed principally for two rea- 
sons. First, to find the effect of the metachronal waves on the flow created and second, to 
find the physical origin of the metachronal waves. Theoretical and numerical studies have 
been undertaken by biologists and fluid m echanicians to understand the flow created by 



an array of cilia (s ee for e.g. the reviews by 



Brennen fe Winet 



1977 



Smith et al 



Blake fc Sleigh 



1974; 



2008). Most of these analyses have been performed to model the flow of 
specific biological systems (e.g. micro-organisms or airway cilia), however, a systematic 
study is lacking. In the following, we outline a number of studies in which the effect 
of the metachronal waves on fluid transport has been studied. Modelling approaches to 



understand t he cilia-driven flow inc lude the envelope model ( Brenne n fc Winet 



Blake 



Liron 



1971a 



1978 



b ), the sublayer model ( Blake 



Gauger et al 



2009 



1972 



Gueron et al 



Gueron fc Levit-Gurevicb 



tion models using a lattice-Boltzmann a 



boundary method (jDauptain et al 



pproach 



Kim fc Netz 



1997 



Smith et al 



1977 



2007 



1999T). flu id structure interac- 



20061 ). and the immersed 



2008f ). In the envelope model, the cilia are assumed 
to be very densely spaced so that the fluid experiences an oscillating surface consisting 
of the tips of the cilia. The envelope model is accurate only when the cilia are spaced 
very close toget her, which has only been observed in the case of symplectic mctachrony 



(Blake 



1971a 



g). In the sublayer model (jBlake 



19721 ). the cilia are represented by a dis- 



tribution of Stokeslets with appropriate mirror images to satisfy the no-slip condition 
on the surface to which the cilia are attached. The sublayer model predicts that for an 
organism that exhibits antiplectic metachrony, the flow created is lower than for cilia 
beating in-phase. In the case of an organism exhi biting symplec t ic me tachrony, the op- 



posite trend is observed. In the numerical study of 



Gauger et al 



( 20091 ). the flow due to 
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the out-of-phase motion of a finite number of magnetic cilia subjected to an oscillating 

external magnetic field was studied. The magnetic cilia generate an asymmetric motion 

due to the difference in the speed of oscilla t ion of the magnetic field during the effective 



and recovery strokes. In contrast to 



Blakd (|l972l) . it was predicted that the flow in the 



case of antiplectic metachrony is larger than the flow created by a symplectic metachrony 
for a particular inter-cilia spacing. 

Early experiments indicated that the hydrodynamic coupling between c ilia could be the 
cause for the formation of the metachronal waves (see for e.g. the review by|] 



Kinosita fc Murakami 



19671) . By mimicking the ciliary motion of Paramecia using an internal actuation mecha- 



nism, it was demonstrated that cilia, which were initially beating in-phase, will form an 



19971) . This behaviour 



antiplectic metachronal wave after a few beat cycles (jGueron et al. 
was explained to be an outcome of the hydrodynamic interactions between neighbouring 
cilia. Similar hydrodyna mically-caused met achronal motion of the cilia was al so observed 
in the numerical work of 



Mitranl (|2007l ). In 



Gueron fc Levit-Gurevichl ((1999), it was re- 



ported that in the presence of the metachronal wave the cilia become more efficient in 
creating flow. The synchroni zation and phase locking of the cili a have also been analysed 



using simple experimen tal (jQian et al 



Vilfan fc Jiilicher 



20091 ) and analytical (jNiedermaver et al 



2008 



20061 ) models. It was found that some degree of flexibility is required 



for t he phase locking of the cilia to take place (INiedermaver et 



2008: 



Qian et al. 



20091 ). The requiremen t of the flexibility fo r synchronization is also confirmed from the 



more detailed model of iKim fc Netd (|2006l ). In the aforementioned studies, however, the 
metachronal wave is an outcome of that specific system, and the flow or the efficiency 
has not been studied for different types of metachronal waves. 

The goal of this paper is, therefore, to obtain a full understanding of the dependence 
of flow on the magnetically-induced out-of-phase motion of an array of asymmetrically 
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beating plate-like artificial cilia at low Reynolds numbers. We will answer the following 

questions using a coupled solid-fluid magneto-mechanical computational model. How 
does the generated flow in the presence of metachrony differ from the flow generated by 
cilia that beat in-phase? How does the flow depend on the metachronal wave speed and 
its direction, and how does it depend on the cilia spacing? We answer these questions 
in the light of magnetic artificial cilia which exhibit an asymmetric motion and beat 
out-of-phase. However, the results are equally applicable to any ciliary system in which 
the cilia exhibit an asymmetric and out-of-phase motion. 

The paper is organised as follows. The boundary value problem, the governing equa- 
tions and the numerical solution methodology are explained in section[2] In section^ the 
physical mechanisms responsible for the enhanced flow in the presence of metachronal 
waves are discussed. The quantitative variation of the flow as a function of the phase 
difference and cilia spacing is given. Finally, the outcome of the analysis is summarised 
in section [U 

2. Problem statement and approach 

We study the flow in an infinitely long channel of height H created by a two-dimensional 
array of plate- like magnetic artificial cilia (having length L and thickness h), which are 
actuated using a rotating magnetic field which is uniform over each cilium, but with a 
phase difference between adjacent cilia. The external magnetic field experienced by the 
i th cilium is 

B x i = B cos(wt - Byi = B a sin(wt - fc), (2.1) 

where Bq is the magnitude of the applied magnetic field, the phase of the magnetic field 
4>i = 2ir(i — l)/n, L) = 2 / n/t- [C i is the angular frequency and i ro f is the time period of 
rotation of the magnetic field. The magnetic field experienced by the individual cilia 
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(a) (b) 
Figure 1. (Color online) (a) Schematic representation of the problem analysed. We study an 
infinitely long microfluidic channel consisting of equal-sized cilia (having length L and thickness 
h) spaced a distance a apart. The variation of magnetic field in space is shown using blue arrows. 
Q p and Q n denote the flow in the direction of the effective and recovery stroke, respectively, (b) 
Typical asymmetric motion of a cilium. The dashed lines represents the trajectory of the tip of 
an individual cilium. 

during a particular instance in time is shown using the blue arrows in Fig. [TJa). The 
phase difference in the applied magnetic field between adjacent cilia is A(f> = lix jn. The 
chosen form of the phase <f>i makes the phase of the magnetic field at every n th cilium 
identical. That is, the magnetic field is periodic after n repeats of cilia. Consequently, 
the applied magnetic field travels n cilia units in time i re f, so that the phase velocity of 
the magnetic field is n/t le { = w/A0 (in cilia per second). The phase velocity is to the 
right (positive) and the magnetic field at each cilium position rotates counterclockwise 
with time. The typical asymmetric motion of a cilium is shown in Fig. Gib). The cilia are 
tethered at one end to the surface, while the other end is free. The trajectory of the free 
end of a typical cilium is represented by the dashed lines in Fig. [Ub), with the arrows 
representing the direction of motion. 

Due to the super-paramagnetic (SPM) nature of the cilia, for which the magnetization 
is proportional to the magnetic field, the magnetic body couple (N = M x Bo, where M 
is the magnetization of the cilia and Bq = {B Xl B y ) is the magnetic field experienced by 
the cilia) depends only on the orientation and magnitude of the magnetic field, but not on 
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its sign. As a result, the bo dy couple at the i t h cilium N Z i, which determines its motion, 



scales with sin (2u)t — 2<fii) (IRoper et 



20061 ). This has consequences for the motion of 



the cilia, both temporally and spatially. Temporally, the frequency of the magnetic couple 
is twice that of the applied magnetic field. This results in two cilia beats for one 360° 
rotation of the magnetic field. Spatially, the phase of the magnetic couple is twice that 
of the applied magnetic field, so that the phase difference between neighbouring cilia is 
twice as large. This means that the magnetic couple is periodic after n/2 cilia. Since 
both the frequency and phase difference increase by a factor 2, the phase velocity of the 
magnetic torque remains equal to that of the magnetic field, i.e. w/A</>. Note, however, 
that the phase velocity of the magnetic torque is equal to the velocity of the metachronal 
wave (i.e., the actually observed dcformational wave travelling over the cilia) only when 
the phase difference A</> is small (i.e. n is large). 

When the phase difference is too large, the metachronal wave can change sign, so that 
the metachronal wave is observed to travel in a direction opposite to the direction of the 
magnetic field (see appendix [AJ. The metachronal wave velocity is equal to u>/ A(j> (i.e. to 
the right) when < A<f> < 7r/2, and it is equal to —u/(ir — A</>) (i.e. to the left) when 
7r/2 < A<f> < 7r, see Fig. [2j When A(j> = 0, the magnetic couple is uniform and all cilia 
beat in-phase. When A<p = 7r, the magnetic couple acting on two neighboring cilia is 
the same (because the phase difference of the magnetic couple is 2A0 = 2ir), and again, 
all the cilia beat in-phase. When A(j> ~ tt/2, the positive metachronal wave velocity is 
equal in magnitude to its negative counterpart. In such a condition, a standing wave is 
observed which causes the adjacent cilia to move in anti-phase. When < A<fi < tt/2 
the metachronal wave velocity is positive, i.e. to the right in Fig. [T] Consequently, the 
metachronal wave velocity is opposite to the direction of the effective stroke, which 
is commonly addressed as antiplectic metachrony (AM). When tt/2 < A<f> < tt, the 
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Figure 2. Metachronal wave velocity as a function of the phase difference A(j> in the mag- 
netic field between adjacent cilia. AM and SM refer to antiplectic and symplectic metachrony 
respectively. 

metachronal wave velocity is in the same direction as the effective stroke and is referred 
to as symplectic metachrony (SM) , see Fig. [2j 

2.1. Governing equations 

We now briefly discuss the coupled solid-fluid magneto-mechanical numerical model used 
to study fluid propulsion using magnetically actuated plate-like artificial cilia. In typical 
microfluidic channels the height H is smaller than the out-of-plane width. Moreover, the 
artificial cilia under study are plate- like (having an out-of-plane width b much larger than 
their thickness h and length L) exhibit a planar beat motion. Therefore, any variation in 
the out-of-plane direction can be neglected and under these assumptions it is sufficient 
to model the artificial cilia and the resulting flow in a two-dimensional setting. 

2.1.1. Solid dynamic model 

We model the cilia as elastic Eulcr-Bcrnoulli beams taking into consideration geometric 
non-linearity in an updated Lagrangian framework. As a starting point for the Euler - 



Bernoulli beam element formulation we use the principle of virtual work (jMalvern 



19771) 
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and equate the virtual work of the external forces at time t + At {SW^ t At ) to the internal 



work (SW^ At ). The internal virtual work is given by 



(aSe + p(uSu + vSv)) dV, 



(2.2) 



where u and v are the axial and transverse displacements of a point on the beam and p 
is the density of the beam. Furthermore, a is the axial stress and e is the corresponding 
strain, given by 



di 



1 / dv 



e dx 2 lar/ V dx 2 ' 



The external virtual work is 



<M& At = / (fxSu + f y 5v + N z ^-) Adx 

J V ° X J (2.3) 

+ J (t x Su + tySv) bdx, 
where f x and f y are the magnetic body forces in the axial and transverse directions, 
N z is the magnetic body couple in the out-of-planc direction, t x and t y are the surface 
tractions and b is the out-of-plane thickness of the cili a. 
We follow the approach used in 



Annabattula et al 



(120101 ) to linearise and discretise 



the principal of virtual work to get, 

5p T (KAp + Mp t+At - Flt t At + F\ nt ) = 0, (2-4) 

where K is the stiffness matrix that combines both materi a l and geometric contributions, 
M is the mass matrix that can be found in 



Cook et al 



(|200lD . F*+ At is the external 



force vector, -F int is the internal force vector, Ap is the nodal displacement increment 
vector and p is the nodal acceleration vector. The nodal acceleration vector is discretized 
in time using Newmark's algorithm (using parameters 7 = 1.0 and (3 = 0.5) so that 
Eqn. 12.41 can be written in terms of the velocity of the beam. The complete discretized 
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equat ions of motion for the solid mechanics model can be found elsewhere (jKhaderi et al 



2009). 



2.1.2. Magneto statics 

To find the resulting magnetic forces, the magnetization of the cilia has to be calculated 
by solving the Maxwell's equations in the deformed configuration at every time increment. 
The Maxwell's equations for the magnetostatic problem with no external currents are 

V-B = VxH = 0, (2.5) 

with the constitutive relation B = /j,q(M + H), where B is the magnetic flux density 
(or magnetic induction), H is the magnetic field, M is the magnetization, and fio is 
the permeability of vacuum. Equatio n 12.51 is solved for M and B using the boundary 



element method ( Khaderi et al 



120091 ) . The magnetic couple per unit volume is given by 
N = M x Bq. As the simulations are two dimensional, the only non-zero component of 
magnetic body couple is N z which is the source for the external virtual work in Eqn. 12.31 
Since the applied magnetic field is uniform for each cilium, the magnetic body forces due 
to field gradients are absent. 

2.1.3. Fluid dynamics and solid fluid coupling 

We study the flow created by artificial cilia in the limit of low Reynolds number. The 
fluid is assumed to be Newtonian and incompressible. The physical behaviour of the fluid 
is governed by the Stokes equation: 

- Vp + 2/iV • D = 0, 

(2.6) 

V -u = 0, 

where p is the pressure in the fluid, D is the rate of deformation tensor, u is the velocity 
of the fluid and fi is the viscosity of the fluid. The set of equations in Eqn. 12.61 is solved 
using Eulcrian finite elements based on the Galerkin method. The fluid domain is dis- 
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cretized into quadrilaterals in which the velocity and pressure of the fluid are interpolated 

quadratically and linearly, respectively. The velocity is calculated at the vertices, mid- 
sides and mid-point of the quadrilateral, and the pressure is calculated at the vertices. 
The solid and fluid domains are coupled by imposing the constraint that the velocity 
at the nodes of the solid beam are equal to the velocity of the surrounding fluid (point 
collocation method) . This coupling is established with the help of Lagrange multipliers 
using the fictitious domain method. Details of the Eulerian finite element model and the 
coupling procedure can be found in Ivan Loon et al\ (|2006f ). 



The fluid domain used for the simulations has a width W and height H (Fig. [3]). For 
each value of a/L, we choose n to be a fraction p/q larger than 2, with p and q integers, 
yielding a range of phase differences A(f> = 2-njn between and 7r. For each value of 
p/q, a unit-cell of width W = pa needs to be chosen to account for periodicity in the 
magnetic couple, unless p is an even integer, for which W = pa/2 suffices. For example, 
let p = 10 and q = 3. Now, n = 10/3 and the phase difference A(j> is equal to 3tt/5. 
To maintain periodicity in the magnetic couple, the width of the unit-cell should be 5a 
(containing 5 cilia). The top and bottom of the unit-cell are the channel walls, on which 
no-slip boundary conditions are applied, 

^top — ^bottom 0, 

while the left and right ends are periodic in velocity 

"left = Wright' 

2.1.4. Solution procedure 

The solution procedure is as follows. The Maxwell's equations are solved at every 
time instant to solve for the magnetic field. From the magnetic field, the magnetic body 
couple acting on the cilia is calculated and is provided as an external load to the coupled 
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Figure 3. Fluid (black) and solid (red) mesh used for the simulations. The mesh corresponds 

to A<fi = 7r/6 and a — 1.67L. 

solid-fluid model, which simultaneously solves for the cilia velocity, and the velocity and 
pressure of the fluid. The velocity of the cilia is integrated using Newmark's algorithm 
to obtain its new position, and the procedure is repeated. Each cilium is discretized into 
40 elements and every fluid domain of size a x H is discretized into 28 x 30 elements, 
with the mesh being refined near each cilium. A typical mesh used for the simulations 
is shown in Fig. [3] A fixed time-step of 1 /is was used for all the simulations reported 
in this paper. The spatial and temporal convergence of the numerical model is discussed 
in appendix [Cj The particles and streamlines are obtained fr om the velocity field in the 



2008h . Also here care should be 



fluid using the visualization software Tecplot (|Tecplot 
taken to accurately resolve the velocity field. 

2.2. Parameter space 

The physical dimensionlcss numbers that govern the behavior of the system are the 
magneto-elastic number M n = V2B 2 L? / figEh 2 - the ratio of the magnetic to the elastic 
forces, the fluid number F n = Yl[il? / Eh 3 tbeat - the ratio of viscous forces acting on the 
cilia to the elastic forces, and the inertia number /„ = 12pL 4 / Eh 2 t?^ t - the ratio of 



Khaderi et al 



20091 ). Here, E is 



the inertia forces of the cilium to its elastic forces, (see 
the elastic modulus of the cilia, h is the thickness, p is the density of the cilia, /i is the 
fluid viscosity, ibcat(— t IC {/2) is the time period of one beat cycle and is the magnetic 



permeability. The geometric parameters that govern the behavior of the system are the 
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phase difference A<f>, the cilia spacing a, their length L and the height of the channel H. 

We study the flow created as a function of the cilia spacing a (normalised with the length 

L) and the phase difference A<fi for the following set of parameters: F n — 0.15, M n = 12.2, 

/„ = 4.8 x 1CP 3 and H/L = 2. The values of the physical parameters correspond to L — 

100 microns, E = 1 MPa, the thickness of cilia being h = 2 (im at the fixed end and 1 

fim at the free end, p = 1600 kg/m 3 , /i = 1 mPas, Bq = 22.6 mT and the cycle time 

i re f = 20ms. The m agnetic susceptibilit ies of the cilia are 4.6 along the length and 0.8 



along the thickness (jvan Riisewiik 



2006) 



The fluid propelled is characterised by two parameters: the net volume of the fluid 
transported during a ciliary beat cycle and the effectiveness. The horizontal velocity field 
in the fluid at any x position, integrated along the channel height gives the instantaneous 
flux through the channel. This flux integrated in time over the effective and recovery 
stroke gives the positive (Q p ) and negative (Q n ) flow, respectively (see Fig. [T]). Due to 
the asymmetric motion, the positive flow is larger than the negative flow, generating a 
net area flow per cycle (Q p — Q n ) in the direction of the effective stroke. The effectiveness, 
defined as (Q p — Q n )/(Qp + Qn), indicates which part of the totally displaced fluid is 
effectively converted into a net flow. An effectiveness of unity represents a unidirectional 
flow. 



3. Results and discussion 

To obtain an understanding of fluid flow due to the out-of-phase motion of cilia, we 
analyse the case of antiplectic metachrony with a phase difference A<f> = 2ir/n = 2-k/YI. 
Since n is even, a unit-cell of width 6a consisting of 6 cilia is chosen, see Fig. |4] The 
contours represent the absolute velocity normalised with L/ibeat- The direction of the 
velocity field can be determined from the arrows on the streamlines. The white arrows 
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represent the applied magnetic field for each cilium. Animations of the ciliary motion for 

the cases of symplectic, antiplectic and anti-phase motion are provided as supplementary 

information. 

The snapshots shown in Figs. 4(a)||4(f)' correspond to the time instances when the flux 



generated by the cilia is maximum. In Fig. 4(g) the instantaneous flux as a function 
of time t (right axis) in addition to the flow (accumulated flux at time t, left axis) are 
plotted. The time instances corresponding to Figs. 4(a)|4(f) are marked in Fig. 4(g) The 
motion of the fluid particles near the third cilium under the influence of the velocity field 



caused by the ciliary motion is also shown. It can be observed from Fig. 4(g) that one beat 
cycle consists of six sub-beats, which correspond to the traveling of the magnetic couple 
from one cilium to the next. The traveling of the metachronal wave to the right can, for 
instance, be seen by looking at the cilia which exhibit the recovery stroke (i.e. cilium 1 
in Fig. 4(a) cilium 2 in Fig. |4(b)| etc). The negative flow created by the cilia during 
their recovery stroke is overcome by the flow due to the effective stroke of the rest of 
the cilia; this leads to a vortex formation near the cilia exhibiting their recovery stroke. 
As a result, the negative flow is completely obstructed for most of the time during the 



recovery stroke. It can be observed from Fig. 4(g) that no flux (right axis) is transported 
in the negative direction, and that the flow (left axis) continuously increases during 
each sub-beat. Moreover, the increase in the flow during each sub-beat is similar (see 



Fig. 4(g)). Thus, the total flow per beat cycle (left axis of Fig. 4(g)) is the sum of the 
flows generated during each sub-beat (i.e. flow per beat = 6 x flow generated during one 
sub-beat). Therefore, it is sufficient to analyse the fluid flow during one sub-beat. 

In the following, we analyse the fluid motion and the resulting flow during the sec- 
ond sub-beat. The velocity profiles at different instants of this sub-beat are shown in 
Figs. 5(a)|5(dJ The corresponding flow and the flux generated are shown in Fig. 5(e) At 
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(f) * = 5tbeat/6 

Figure 4. (Color online) (a)-(f) Out-of-phase motion of cilia during a representative cycle for 
Acj> = 7r/6 (n — 12) with the wave moving to the right (antiplectic metachrony) for a/L — 1.67. 
The contours represent the absolute velocity normalised with L/ibeat- The direction of the 
velocity is represented by streamlines. The white circles represent fluid particles. The applied 
magnetic field at each cilium is represented by the white arrows, (g) Instantaneous flux (right 
axis) and flow (or accumulated flux, left axis) as a function of time with the instants (a)-(g) 
duly marked. 
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ibeat/6, the third cilium starts its recovery stroke and the particles near the top bound- 



ary are driven by the positive flow created by cilia 4, 5 and 6 (sec Fig. 5(a)). At this 
instant, as only one cilium is exhibiting a recovery stroke, the flux created by the cilia is 



maximum (see instant 'a' in Fig. 5(e) I . In Fig. 5(b) the third cilium also has begun its 
recovery stroke and now the negative flow caused by both the second and third cilia is 
opposed by the effective stroke of the other cilia. The high velocity of the second cilium 
during its recovery stroke decreases the flux caused by the other cilia (see instant 'b' in 



Fig. 5(c) ). When the third cilium is half-way through its recovery stroke (see Fig. 5(c) ), 
the second cilium is a bout to finish its rec overy, which generates a large velocity, due to 



the whip-like action (jKhaderi et al 



20091 ) . to the right. Now, the position of the third 
cilium is such, that it opposes the negative flow caused by the second cilium. This leads 
to a strong vortex formation near the second and third cilia, with only a small flux in 
the direction of the recovery stroke (to the right). The small negative flux caused by the 



whip-like motion of the second cilium can be seen by the instant marked 'c' in Fig. 5(e) 



causing a momentary decrease in the flow. The vortex imparts a high velocity in the 
direction of the effective stroke to the particles away from the cilia. As the third cilium 
progresses further in its recovery stroke, the particles come under the influence of the 
flow due to the rest of cilia, which are now in different phases of their effective stroke 



(see Fig. 5(d)). Now, only the third cilium is in the recovery stroke; this again leads to 



a maximum value of the flux (similar to Fig. 5(a) ). The key observation of Figs. [4] and 
O is that the negative flow created during the recovery stroke of the cilia creates a local 
vortex due to the positive flow created by other cilia. This shielding effect during the 
recovery stroke leads to a drastic increase in the net propulsion rate for cilia beating 
out-of-phase, compared to synchronously beating cilia. 

Next, we analyse the instantaneous flux (Fig. 6(a) ) and flow generated (Fig. |6(b)| as 
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Figure 5. (Color online) (a)-(d) Snapshots for the out-of-phase motion of cilia between time 
instances of Figs. |4(b)| and |4(c)| for = ty/6 (n — 12) with the wave moving to the right 
(antiplectic metachrony) for a/L = 1.67. The contours represent the absolute velocity normalised 
with L/tbe&t- The direction of the velocity is represented by streamlines. The white circles 
represent fluid particles. The applied magnetic field at each cilium is represented by the white 
arrows, (e) Instantaneous flux (right axis) and flow (left axis) as a function of time with the 
instances (a)-(d) duly marked. 

a function of time for different phase differences. When the cilia move synchronously 



(A(j> = 0), the flux (see the solid line in Fig. 6(a)) is positive for approximately three- 



quarters of the time and strongly negative during the rest of the cycle. Consequently, the 
flow generated (see the solid line in Fig. |6(b)[ ) increases during the effective stroke, but 
profoundly decreases when the recovery stroke takes place. This creates a large fluctuation 
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(a) 




(b) 

Figure 6. (a) Normalised fluid flux as a function of time for a/L = 1.67 and different values of 
phase difference A</>. (b) Normalised accumulated flow at any time t during the beat cycle. 

in the flow, with only a small net amount of fluid transported. Once the ciliary motion 
is metachronal, the negative flux is very small compared to the positive flow (see the 



cases of a standing wave and antiplectic metachrony in Fig. 6(a)). This decreases the 
fluctuation in the flow generated, causing it to increase nearly monotonously during the 
beat cycle (see the dashed and dotted lines in Fig. |6(b)| . We can clearly see that the 
flow at the end of the beat cycle (t = ibeat) for out-of-phase motion is significantly larger 
than the flow created by the synchronously beating cilia. 
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The fluid propelled and the corresponding effectiveness are plotted for different values 

of A(j) and a/L in Fig. [7J The mctachronal wave velocity (Fig. [5]) is plotted as a function 



of A</> and is shown using dashed lines in Fig. 7(a) As mentioned earlier, when the 
metachronal wave velocity is positive an antiplectic metachrony (AM) results, and when 
the metachronal wave velocity is negative we get a symplectic metachrony (SM). When 
all the cilia are moving synchronously (A<f> = or it), the flow (normalised by 7rL 2 /2) will 
be approximately 0.22 for a/L = 5. As the cilia density is increased by decreasing a from 
a/L = 5 to a/L = 1.67, the viscous resistance per cilium decreases, which causes the 
normalised flow to increase to 0.25. When the cilia beat in-phase, the effectiveness of fluid 
propulsion is very low, see Fig. |7(b)] The fluid propelled shows a substantial increase once 



the cilia start beating out-of-phase (Fig. 7(a) ). When the cilia spacing is large (a/L = 5 
and 2.5), the flow generated remains approximately constant for all metachronal wave 
speeds. The increase in flow by decreasing the cilia spacing from a/L = 5 to a/L = 2.5 is 
much larger when the cilia beat out-of-phase compared to the increase when the cilia beat 
in-phasc. However, when the cilia spacing is low (a/L = 1.67), we see a larger increase 
in the fluid flow when there is an antiplectic metachrony(AM) compared to a symplectic 
metachrony (SM). Also, the effectiveness sharply increases from around 0.3 (i.e., 30% 
of the totally displaced fluid is converted into net flow) to 1 (fully unidirectional flow), 
sec Fig. |7(b)| To analyse these trends a bit further, we plot the positive and negative 
flow (Q p and Q n in Fig. [I]) created during a beat cycle for different phase differences 



in Fig. 8(a) It can be seen that the cilia do not create a negative flow when they beat 
out-of-phase for all cilia spacings, resulting in a unidirectional flow (effectiveness = 1). 
This reduction in negative flow is due to the shielding of flow during the recovery stroke 
caused by the effective flow of other cilia. It can also be noted that the positive flow is 
also reduced compared to in-phasc beating, but the reduction is considerably less than 
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Figure 7. Flow and effectiveness as a function of the phase difference A<f> for different inter-cil- 
ium spacings a/L. AM and SM refer to antiplectic metachrony (the wave direction is opposite 
to the direction of the effective stroke) and symplectic metachrony (the wave direction and the 
effective stroke direction are the same), respectively. 

the reduction in negative flow. Thus, the net flow increases as soon as the cilia start 



to beat out-of-phase (see Fig. 7(a) I. It can be seen from Fig. 8(a) that in the presence 
of metachronal waves when the cilia spacing is large (a/L = 5), the fluid transported 
during the effective stroke remains nearly the same for all values of the wave velocities. 
For small cilia spacing (a/L = 1.67), however, the positive flow is maximal for antiplectic 
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metachrony, which leads to a larger net flow for antiplcctic metachrony compared to 

symplcctic metachrony. 

To understand the difference in positive flow for opposite wave directions for small 

inter-cilium spacing (a/L = 1.67), we plot the flux as a function of time scaled with 

the time taken by the magnetic couple to travel from one cilium to the next t\, for two 

different metachronal wave velocities (3/£beat and 6/ibeat cilia per second), see Fig. |8(b)| 

The corresponding phase differences are also shown in the legend. It can be seen that 

the flux in the case of antiplcctic metachrony is larger than the flux created by the 

symplectic metachrony for the same wave speed. This difference in flux for opposite wave 

directions can be understood by analysing the velocity field corresponding to symplectic 

and antiplectic metachrony at time instances when the flux is maximum (see Fig. [9]). 

FigureEJa) and[9jb) correspond to different phase differences (Ac/) = ir/6 and A<fi = 5ir/6, 

respectively) leading to a similar wave speed of 6/ibeat cilia per second (see also Fig. [5]). 

The fifth cilium is in the peak of its effective stroke for both AM and SM. In the case 

of symplectic metachrony, the positive flow created by the fifth cilium is obstructed by 

the close proximity of the fourth cilium, which has just started its effective stroke. As 

a result, we observe the formation of a vortex. In the case of antiplectic metachrony, 

however, the position of the fourth cilium is such that the positive flow created by the 

fifth cilium is not obstructed. This leads to larger fluid flow in the positive direction, so 

that the net flow created by an antiplcctic metachrony is larger than that created by its 

symplectic counterpart. 

Reports on metachrony and phase locking of beating cilia have appeared in the past 



( Gaueer et al 



2009 



Kim fc Netz 



19991 



2006: 



Gueron et al. 



1997 



Gueron fc Levit-Gurevicb 



19991 ). The main results are that metachrony enha nces flow compared to synchronously 



beating cilia (|Kim fc Netz 



2006 



Gauger et al 



20091 ) and that antiplectic metachrony gen- 
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Figure 8. (a) Positive (Q p ) and negative flow (Q n ) (see Fig. [1} created by the cilia corre- 



sponding to the results presented in Fig. [7] (b) Flux vs time (scaled with the time t\ taken by 
the magnetic couple to travel from one cilium to the next) for a/L = 1.67 and different wave 
speeds. 



erates a higher flow rate than symplectic metachrony (jGauger et al 



2009) 



Kim&Netz 



( 2006T ) analysed two cilia, which are driven by internal motors and are moving out-of- 
phase due to the hydrodynamic interaction. They have shown that the fluid propulsion 
increases, once the cilia start to beat wi th a phase difference , which is in agreement with 



our results. Our results also agree with 



Gauger et al 



(2009), where it is shown that the 
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1 2 3 A 5 6 

(a) Antiplectic metachrony: wave travels to the right 




(b) Symplectic metachrony: wave travels to the left 
Figure 9. (Colour online) Snapshots for antiplectic (A(fi = 7r/6) and symplectic metachrony 
(A<f) — 5n/6) for a wave speed of 6/tbcat cilia per second and cilia spacing a/L = 1.67 at t — O.lti 



of Fig. 8(b) The contours represent the absolute velocity normalised with L/tbeat (blue and red 



colours represent a normalised velocity of and 2, respectively). The direction of the velocity is 
represented by streamlines. The applied magnetic field is shown by the white arrows. 

fluid flow is larger in the case of antiplectic metachrony than s ymplectic metachron y when 



Gauger et al 



(2003) in the 



the cilia are close together. However, our results differ from 
sense that we always see an enhancement in flow in the presence of metachrony (compared 
to cilia beating in-phase) irrespective of the direction and magnitude of the metachronal 
wave velocity. This is most likely due to the fact that the asymmetry in ciliary motion 



in our case is much higher. 



Gueron et al 



( 1997 ) and 



Gueron fc Levit-Gurevichl ( 19991) 



have proposed that the evolution of the out-of-phase motion of cilia in Paramecia is due 
to hydrodynamic interactions between adjacent cilia leading to antiplectic metachrony. 
It is interesting to observe that the interplay between the internally-driven actuation 
and hydrodynamic interaction in nature results in antiplectic metachrony Our results, 



and those of others ( Gauger et al 



120091 ) , show that indeed antiplectic metachrony leads 
to larger flow than symplectic metachrony for small cilia spacings as typically seen in 
nature. 
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4. Conclusions 

Using a numerical model we have studied the flow created by a two-dimensional array 
of plate-like artificial cilia as a function of the phase lag and spacing between neighbouring 
cilia. The flow per cycle and the effectiveness (which is a measure of the unidirectionality 
of flow) are considerably enhanced when the cilia start beating out-of-phase, as compared 
to synchronously beating cilia. While the amount of flow enhancement depends on the 
inter-cilia spacing, the effectiveness is not significantly influenced. Metachrony is observed 
to completely knock-down the negative flow to zero due to the vortex formation caused 
by the shielding of the recovery stroke. Interestingly, we find that the enhancement is 
achieved even for small phase differences. The direction of travel of the metachronal wave 
is important only for small cilia spacing. In that case, the flow is larger for antiplectic 
metachrony compared to symplectic metachrony, which is related to the obstruction 
of the positive flow for symplectic metachrony. It is therefore beneficial if the magnetic 
actuation of the artificial cilia is designed such that it results in an antiplectic metachrony. 
Our results suggest that an antiplectic metachrony is adopted by the cilia on paramccia 
and in the respiratory system to maximize the fluid propelled. However, ciliary systems 
(such as on Opalina) that exhibit symplectic metachrony are also present in nature. It 
will be of interest to investigate what property is optimised by symplectic metachrony in 
these systems. 
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Appendix A. Metachronal wave velocity 

The metachronal wave velocity is obtained by dividing the distance between two cilia 
with the time it takes for the magnetic couple to travel from a cilium to its neighbor. If 
the neighbor is to the right, then the wave travels to the right, and when the neighbor 
is to the left, the wave travels to the left. The magnetic couple Ni at any cilium i is 
proportional to sin (2ujt — 2<f>i), and travels with a phase velocity of w/A</> (in number of 
cilia per second) to the right. 

In the schematic of Fig. [lUJ three cilia C\. C2 and C3 are depicted. At any given 
instance of time, let the magnitude of the magnetic couple at Ci, C2 and C3 be Ni, N2 
and N3, respectively. The magnitude of the magnetic couple at the 'periodic' cilium H, 
which is separated from C3 by n/2 units, is also N3. The metachronal wave is said to 
have traveled to the right when the magnetic field at C2 is Ni after a time interval. Now, 
the distance traveled by the magnetic couple is 1 cilia spacing, and the time taken to 
travel this distance is l/(w/A0). Therefore, the velocity of the magnetic couple is w/A</>, 
in cilia units per second. The metachronal wave is said to have traveled to the left when 
the magnetic field at C2 is equal to _/V 3 after an interval of time. As the applied magnetic 
couple travels to the right, this situation is possible when the magnetic couple at the 
periodic cilium H travels to the cilium C2. The time needed for the magnetic couple 
to travel from H to C2 is equal to (n/2 — l)/(w/A</>). However, the apparent distance 
travelled is one cilium spacing to the left (i.e. from C3 to C2), so that the wave velocity 
is now lo/(tt — A</>). The (apparent) metachronal wave velocity is now determined by the 
maximum of the two competing wave velocities: w/A0 to the right and u)/(tt — A</>) to 
the left. As a result, the metachronal wave velocity is equal to w/A</> (i.e. to the right) 
when u)/A(j) > oj/(tt — A(j)) (i.e. < A0 < tt/2), and it is equal to —u/(ir — A<f>) (i.e. to 
the left) when w/A(j) < uj/(jz - A0) (i.e. tt/2 < A(f>< tt), sec Fig.0 



Metachronal beating of magnetic artificial cilia 27 



H C, 
nil cilia units 



Figure 10. Schematic diagram used to calculate the metachronal wave velocity. 
Appendix B. Validation of the fluid-structure interaction model 

To compare the performance of the present approach with a solution available in the 
literature we choose to study the deformation behavior of a cantilever beam under an 
imposed pulsating flow. This problem has been numerically solved by iBaaiiena (|2001l) 



using the fictitious domain method in which the solid was discretized using continuum 
finite elements. The width W is four times the height H of the fluid domain. H is taken 
to be unity. The length of the cilium is 0.8H. The thickness of the cilium is 0.0212iJ. The 
elastic modulus of the cilium and viscosity of the fluid were specified in dimensionless 
units to be E = 10 7 and fi = 10, respectively. The mesh used for the computation is 
shown in Fig. [TT] The dots represent the nodes of the Euler-Bernoulli beam element. 
The boundary conditions are as follows: the left and right boundaries are periodic. A 
pulsating flow of magnitude 10 sin(27rf/T) is prescribed on the left boundary, where T 
is the time period which is taken to be sufficiently large to avoid inertia effects in the 
cilium. The bottom boundary is a no slip boundary. On the top boundary, the normal 
flow is constra i ned. T he solution from our formulation is plotted along with the solution 
from 



Baaijensl (|2001l ) in Fig. IT27 a) in terms of the displacement of the free end of the 
cantilever. It can been seen that the two solutions are in good agreement. In Fig. [T27b). 
we plot the x displacement of the free end of the beam as a function of time for different 
discretizations of the cilium (using 12, 24 and 48 beam elements). When the cilium mesh 
is refined, the fluid mesh is also refined proportionally, see also appendix C. It can be 
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Figure 11. Coarsest mesh used for benchmarking. 
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Figure 12. A cantile ver subjec t ed to a pulsating flow: Comparison of solution obtained from 
the present work with 



Baaiiem 



(2001). (a) Comparison of the trajectory of the free end. The 



deformed and initial configurations are also shown, (b) Comparison of the displacement of the 
free end as a function of time for various mesh refinements. The number in parenthesis of the 
legend refer to the number of elements used to discretize the cantilever. 

seen that the displacements nicely converge as the mesh is refined. The convergence of 
the velocity field is also shown in Fig. Q21 



Appendix C. Convergence of the numerical model 

In this section, we report on the spatial and temporal convergence of the numerical 
method used in this paper. We use the case of synchronously beating cilia (A</> = 0) with 
an inter-cilia spacing of a = 1.67L, for which the unit-cell consists of one cilium. As the 
deformed shape of the cilium is an outcome of the model, wc compare the position of 
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Figure 13. Convergence of velocity field at a particular time instant with mesh refinement. The 
mesh used for Fig. (a) is shown in Fig. 1111 where 12 beam elements are used. In Fig. (b) and 
(c), 24 and 48 elements were used to discretize the cilia, while the fluid mesh was also refined 
proportionally. 

the free end for different temporal discretizations. The mesh used to discretize the cilium 
and the fluid domain is shown in Fig. [14] for the case when the cilium is divided into 40 
cilia elements and the fluid is divided into 28 x 30 elements. 

The position of the tip of the cilium as a function of time and its trajectory for different 
time increments is shown in Fig. [15] (a)-(c). The time increment has to be small enough 
to capture the fast whip-like recovery stroke. It can be seen that a time increment of 
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Figure 14. Discretization used for cilium and fluid. The cilium is discretized into 40 elements 
and the fluid domain of size a x H is divided into 28 x 30 elements. 

1 fj,s is sufficient for temporal convergence. This time step of 1 /its is used to study the 

spatial convergence and the results are shown in Fig. [16] The number of elements on the 

cilium as well as the fluid are changed proportionally when the mesh is changed. In the 

following the spatial discretization is defined in terms of the number of elements used to 

discretise the cilium; i.e., 30 cilia elements correspond to a fluid mesh of 21 x 23 and 60 

cilia elements correspond to a fluid mesh of 42 x 45. It can be seen that the results for 

these discretizations have fully converged as shown for the position of the free end of the 

cilium and the flux function of time. 
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